Jamming, hysteresis and oscillation in scalar models for shear thickening 
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We investigate shear thickening and jamming within the framework of a family of spatially 
homogeneous, scalar rheological models. These are based on the 'soft glassy rheology' model of 
Sollich et al. [Phys. Rev. Lett. 78, 2020 (1997)], but with an effective temperature x that is a 
decreasing function of either the global stress a or the local strain I. For appropiate x = x(a), it 
is shown that the flow curves include a region of negative slope, around which the stress exhibits 
hysteresis under a cyclically varying imposed strain rate j. A subclass of these x(a) have flow curves 
that touch the 7 = axis for a finite range of stresses; imposing a stress from this range jams the 
system, in the sense that the strain 7 creeps only logarithmically with time t, y(t) ~ Int. These same 
systems may produce a finite asymptotic yield stress under an imposed strain, in a manner that 
depends on the entire stress history of the sample, a phenomenon we refer to as history-dependent 
jamming. In contrast, when x = x(l) the flow curves are always monotonic, but we show that 
some x(l) generate an oscillatory strain response for a range of steady imposed stresses. Similar 
spontaneous oscillations are observed in a simplified model with fewer degrees of freedom. We 
discuss this result in relation to the temporal instabilities observed in rheological experiments and 
stick-slip behaviour found in other contexts, and comment on the possible relationship with 'delay 
differential equations' that are known to produce oscillations and chaos. 

PACS numbers: 83.60.Rs, 64.70.Pf, 83.10.Gr 



I. INTRODUCTION 

A wide variety of materials can be driven into a non- 
equilibrium state that is either solid-like and static, or 
fluid-like but only relaxes on time scales that far ex- 
ceed the experimental time frame, if at all. Such states 
are often referred to as 'jammed,' and are realisable in 
molecular liquids having undergone a rapid quench to 
low temperatures, or colloids at a high volume fraction, 
to cite just two examples It has recently been pos- 

tulated by Liu and Nagel that the nature of the jammed 
state may be independent of the manner in which it 
was formed Q. For example, a stress-induced jamming 
transition may produce a qualitatively similar state to a 
temperature-induced transition. This was expressed in 
the form of a 'jamming phase diagram,' in which jammed 
configurations occupy a compact region near to the origin 
of a phase space comprised of three axes : the tempera- 
ture T, the volume V and the load a 

In its simplest form, the Liu-Nagel jamming phase di- 
agram suggests that increasing the applied load can only 
increase the likelihood of flow. However, it is equally 
feasible for an applied load to induce a jammed state. 
For instance, if a pile of sand is formed and then gravity 
is switched off, the system unjams without any signifi- 
cant variation in volume. Thus the the load (here con- 
trolled by gravity) jams the system. This concurs with 
the earlier claim that jammed systems may be classified 
as 'fragile' — that is, they can support only certain, com- 
patible loads, and will rearrange or flow under an incom- 
patible load g. Furthermore, one could also envisage a 
class of loadings that can alternately induce and destroy 



a jammed state as the magnitude of the load increases, 
a situation that could be referred to as 're-entrant jam- 
ming.' 

The purpose of this paper is to investigate the na- 
ture of transitions to or from jammed configurations 
that have, as their control parameter, the imposed shear 
stress a. We do this by providing concrete examples 
of models that exhibit jamming transitions with er, in- 
cluding instances of the re-entrant jamming scenario de- 
scribed above. These models are based on the 'soft glassy 
rheology' (SGR) model of Sollich et al. ||-^), which was 
originally devised to highlight the possible existence of 
glassy relaxation in a range of soft materials, such as 
foams, emulsions, pastes etc. It is parameterised by an 
effective temperature x, which in in the context of soft 
matter is thought to represent some form of mechanical 
noise, but may refer to the true thermal temperature in 
other applications. 

As it was originally defined, the SGR model only ex- 
hibits shear thinning, which is clearly unsuitable for our 
purposes. Therefore we consider variants of the model in 
which the effective temperature x is no longer constant, 
but can vary with the state of the system. More precisely, 
x is treated a function of both the global stress a and 
the local strain I, i.e. x = x(cr,l). By choosing suitable 
forms of x(a, I), systems can be constructed that become 
'colder' as they become more stressed, which allows for 
the system to shear thicken and even 'jam.' 

For clarity, we restrict our attention to two limiting 
forms of x(a, I), namely x{a) and x(l). In the former case, 
certain classes of x{a) are shown to exhibit a flow curve 
[i.e. the curve of the stress a versus the strain rate 7 un- 
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der conditions of steady flow) that is non-monotonic, and 
can touch the 7 = axis for range of finite values of a. 
Thus a 'jammed' state with 7 = + can be reached for 
some a but not for others. Furthermore, it is also shown 
that, for systems driven by an imposed j(t), whether or 
not a jammed configuration is approached at late times 
depends on the entire strain history of the system, and 
not just the behaviour of j(t) as t — > 00. We refer to this 
phenomenon as 'history-dependent jamming.' 

The behaviour of systems with x — x(l) is somewhat 
different, but no less remarkable. When driven by an im- 
posed a, certain forms of x(l) exhibit a regime in which 
the viscosity never reaches its steady flow value but os- 
cillates in time, with a waveform that is approximately 
sinusoidal near to the transition point, but becomes in- 
creasingly sharp deep into the oscillatory regime. The 
models considered here are spatially homogeneous, and 
so this oscillation is purely temporal, having no spatial 
component. It is tempting to associate these oscillations 
with the stick-slip behaviour observed in granular sys- 
tems and plate tectonics ^-Q, but we shall give rea- 
sons why we believe that the underlying physics might 
be rather different. We cannot rule out the possibil- 
ity of more complex oscillatory behaviour, maybe even 
chaotic, arising in as-of-yet unobserved regions of pa- 
rameter space. 

The finding of a bifurcation to oscillatory behaviour for 
x = x(l) is all the more remarkable because the flow curve 
(as already defined) is everywhere monotonic increasing 
for this class of x. By contrast, instances of rheological 
instabilities that have been found experimentally have 
tended to occur for ranges of parameters in which the 
flow curve has a negative gradient fl5|-|l9[] . This suggests 
that mechanism behind the instability observed here may 
be qualitatively different to those cited above, and may 
be realisable in some range of materials that has yet to 
be identified. 

This paper is arranged as follows. In Sec. || the class 
of models to be studied is defined, with particular atten- 
tion being paid to those aspects that differ from the SGR 
model. The known results of the SGR model that are rel- 
evant to the subsequent discussion are then briefly sum- 
marised in Sec. III. Systems with x — x(a) are described 
in Sec. IV, where it explained how the flow curves can be 



graphically interpreted as mappings from the SGR flow 
curves. The time-dependent behaviour of the system 
has been found by numerical integration of the govern- 
ing master equation. An example of history-dependent 
jamming is presented, and explained in terms of the sta- 
bility of the steady flow solutions. 

In Sec. [v]we turn to consider the case x — x(l), and an- 
alytically prove that the flow curve is everywhere mono- 
tonic. Nonetheless, the simulations results presented here 
clearly show that j(t) can oscillate in time for a range 
of imposed stresses a. A qualitative explanation of the 
emergence of the oscillatory phase is also given, in terms 
of the model's internal degrees of freedom. The intensive 
nature of the simulations has meant that only a small 



range of functional forms for x(l) have been investigated, 
and hence it has not been possible to fully characterise 
this range of models. Therefore we consider a simplified 
model in Sec. for which the simulation times are sig- 
nificantly shorter and a more complete picture has been 
realised. Some results of this model have also been pre- 
sented elsewhere . This reduced model clearly shows 
that the mean value of 7 during the oscillatory motion 
deviates from the steady flow value by a s much as an or- 
der of magnitude. Finally, in Sec. VII we discuss some 
outstanding issues rais ed by this work, before summaris- 
ing our results in Sec. VIII. 



II. MODELS OF THE SGR— TYPE 

All of the models studied in this paper represent gen- 
eralisations of the SGR model of Sollich et al. |^-^|. It 
is therefore prudent to first describe those features that 
are common to this class of models, before considering 
each of the different generalisations in turn. This is the 
purpose of the current section. Since the various assump- 
tions that lie behind the construction of the SGR model 
have already been discussed at length elsewhere, we shall 
here give just a brief overview of the derivation and refer 
the reader to for more detailed physical arguments. 
Only those aspects that differ from the SGR model will 
be discussed in full. 

Our goal is to construct deliberately simplified models 
that exhibit shear-thickening and jamming, but whose 
interpretation is more transparent than that of a detailed 
microscopic model. In this spirit, the models in this class 
all share a number of simplifications. Only a single shear 
component of the strain and stress tensors are considered, 
which will be denoted by 7 and a respectively. This is 
therefore a class of scalar models. 

It is further assumed that the system can be coarse 
grained into a collection of mesoscopic subsystems, each 
of which are fully described by two scalar variables, 
namely a local strain I and a stability parameter E. 
These mesoscopic subsystems will be referred to as el- 
ements. We simply assume here that such a coarse- 
graining is possible, notwithstanding the significant prac- 
tical challenges in constructing a suitable scheme for any 
given microscopic model. 

At any given instant, each element has a probability 
of yielding per unit time that is denoted by T = T(E, I). 
When an element yields, it is assumed that its micro- 
scopic constituents are rearranged to such a degree that 
it loses all memory of its former configuration. Its strain 
returns to zero, and it is assigned a new value of E that is 
drawn from a prior distribution p(E) . Suitable functional 
forms for T(E, I) and p(E) will be discussed below. The 
value of E remains fixed until the element yields again, 
but I follows the global strain 7 according to I = 7. Thus 
the elements are affinely deformed by the flow field, which 
is assumed to be homogeneous. 
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Let us define P(E,l,t) to be the probability density 
function of elements which have a local strain I and a 
barrier E at time t. Then P(E,l,t) evolves in time ac- 
cording to two distinct mechanisms : the homogeneous 
shearing at a rate 7, and the yielding of elements at a 
rate of T(E, I) per unit time. Thus the master equation 
for P(E,l,i) is 



-p(e, i, t) + %W M) = - r(s, op(s, 1, t) 

+ oj(t) 6(1) p{E) . 



(1) 



The second term on the left-hand side represents the in- 
crease in local strains I according to the uniform global 
strain rate 7. The right hand side describes the yielding 
of elements. The first term, which is negative, accounts 
for the loss of elements as they yield at a rate T(E,l). 
Conversely, the second term represents these same ele- 
ments after they have yielded, which have a strain I = 
and a value of E drawn from p(E). The total rate of 
yielding u>(t) is defined by 

poo poo 

w(t)= / dE dlT(E,l)P(E,l,t) (2) 

JO J -00 

= <IW)> • (3) 

Here we have introduced the notation that the angled 
brackets '(■■■)' represents the instantaneous average of 
the given function over P(E,l,t). This will be used fre- 
quently in what follows. 

The master equation (Q) only describes the evolution 
of the strain degrees of freedom. To characterise the rhe- 
ological response of the system, some relation must be 
found between the local strains I and the global stress a. 
This involves two further assumptions. Firstly, the ele- 
ments are supposed to behave elastically between yield 
events, so that the local stress is just kl, where k > is 
a uniform elastic constant. Secondly, a is the arithmetic 
mean of the local stresses, or a = (kl) = k(l). Other av- 
eraging procedures could be employed , but we focus 
here on the simplest non-trivial option. Although the lo- 
cal stress-strain relationship is elastic, the global stress 
also incorporates the yielding of elements and, as will be 
seen below, the cr-7 relationship is not a linear one. 

We have been unable to find an analytic solution to the 
master equation ([l]) for any non-trivial T(E,l). Instead 
it has been numerically integrated using the procedure 
summarised in Appendix |A|. However, it will sometimes 
be necessary to refer to the steady state solution for a 
constant 7 ^ 0, which can be found exactly, 



Poo(EJ) 



7 



-p(E) exp 



1 f l 

— / dl'T(EJ') 
7 Jo 



(4) 



This is derived by setting dtP — Q in (Q) and in- 
tegrating the resulting first-order ordinary differential 
equation with respect to I. The asymptotic yield rate 
LUac = linit_,oo w(t) can be found by requiring that the 
integral of Poo(E, I) is unity. 



A. Yielding modelled as an activated process 



In the SGR model, the yield rate T(E, I) was assigned 
a functional form similar to that of an activated pro- 
cess 0j. This was based on the idea that each element 
can be represented by a single particle moving on a free 
energy landscape, which remains confined within a well 
of depth E until random fluctuations allow it to cross 
over this barrier into a different well, with a new bar- 
rier E' . Thus yielding can be identified with barrier cross- 
ing. This description is similar to that of activated pro- 
cesses in thermodynamic systems; however, the random 
fluctuations here are not necessarily due to the true ther- 
modynamic temperature. They may arise rather from a 
form of homogeneous mechanical noise generated by non- 
linear couplings between the elements; see (Q for a fuller 
discussion on this point. To avoid possible confusion, this 
effective temperature is denoted by the symbol x rather 
than T. 

Although the energy barrier of an element is ini- 
tially E, as it becomes strained it will gain an elastic 
energy of \kl 2 and thus will have a smaller effective en- 



ergy barrier AE 
take the form 



E — ±kl A . Thus the yield rate will 



T(E, I) = To exp 



E 



\kl 2 



(5) 



where the attempt rate Tq sets the time scale of the 
yielding. The effective temperature x is constant in 
the SGR model, and essentially acts as a parameter of 
the model. However, since x may be generated in part 
by internal couplings between the elements, it should 
be allowed to vary with the state of the system, i.e. 
x = x({P(E,l,t)}). 

In this paper, we shall not consider the most general 
possible form for x, but shall instead focus on a more 
restricted class for which x — x(o~,l). This corresponds 
to the realisation that, as an element is strained, it may 
become more or less susceptible to the noise, and hence 
its 'temperature' x may change. Allowing x to also de- 
pend on a reflects that this change in susceptibility to 
noise may in part be a global phenomenon, i.e. the x of 
a given element may depend on the state of all of the el- 
ements around it. Clearly, deriving the actual x(a, I) for 
any given material would be a highly complicated task, 
of comparable difficulty to the original coarse-graining 
procedure described above. 

For further simplicity, we shall not consider x = x(a, I) 
but shall instead focus on two limiting cases : x = x(a) is 
described in Sec. IV, and x = x(l) is assumed in Sees. |v| 
and VI. The relevant results for the SGR model, which 



corresponds to the case of x = (constant), are also sum- 
marised in Sec. III. 
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B. Natural units and the p(E) distribution 



III. SUMMARY OF THE STANDARD SGR 
MODEL 



We have yet to specify p(E), the distribution of en- 
ergy barriers E for elements that have just yielded. In 
the SGR model, p(E) is assumed to have an exponential 
tail, p(E) ~ e~ E / E ° as E — ► oo, which gives rise to a fi- 
nite yield stress and diverging viscosity for some values of 
x, but not others Qj. Although it is possible to justify the 
occurrence of an exponential tail in some contexts [p2[ , 
we prefer instead to treat this choice of p(E), combined 
with the Arrhenius form of T(E,l) (|5|), as a recipe for 
generating a yield stress within this simple picture of ac- 
tivated yielding. This may seem ad hoc, but it should be 
realised that jamming is in reality a many body effect, 
involving collective behaviour between a large number of 
degrees of freedom. It should therefore come as no sur- 
prise that assumptions are required to describe jamming 
within this single-particle picture. 

For much of this paper, we shall assume that p(E) has 
an exponential tail, just as in the SGR model. In fact, 
for the numerics the definite form p(E) — -j^ c ~ E ^ E ° na s 
been used, although it is not expected that the precise 
shape of p(E) for small E will alter the long-time or 
steady state behaviour of the system. This is because 
small values of E correspond to elements with short ex- 
pected lifetimes. Furthermore, it turns out that the os- 
cillatory motion described in Sec. M does not depend 
on the choice of p(E), and in Sec. ^ the simpler form 
p(E) = 8(E — Ei) is used, with qualitatively similar re- 
sults. 

There remain three constants in the problem. These 
are the elastic constant k, the attempt rate Tq in ([|), and 
the constant Eq in the definition of p(E) given above. 
However, these can all be scaled out of the model. For 
instance, k sets the scale of the stress, and therefore can 
be removed by rescaling a to a/k. Similarly, Tq gives the 
only intrinsic time scale in the system, and can be scaled 
out through t and 7. Finally, Eq sets the scale of the 
energy barriers, and can be scaled out through x(a,l). 

In what follows, we have adopted natural units in which 
k = To = Eq = 1. This means that the actual values for 
t, 7 and a given below should not be compared to experi- 
mental values without first rescaling with the appropriate 
k, Eq and Tq for the material in question. For example, 
the typical yield strain according to the T(E, I) given in 
(|) is y / 2E /k, which is when the effective energy barrier 
vanishes. In natural units, this is of 0(1); however, for 
soft materials it will typically be of only a few percent, 
and will be even smaller for hard materials. Thus the 
scale of the local strains I, and therefore the scale of the 
global strain "f(t), will generally be significantly smaller 
in real materials than with natural units. 



As discussed in the previous section, the standard SGR 
model is realised when (using natural units) the yield 
rate T(E, I) is chosen to take an Arrhenius form with an 
energy barrier E — \l 2 and a constant effective temper- 
ature x, and p(E) has an exponential tail p(E) ~ e~ E . 
Many results for this case are already known [| . The 
purpose of this section is to briefly describe those results 
that will be referred to in later sections. 

Let us assume that the system is driven by a given 
strain j(t), and that 7(f) ~ jt as t — > 00, where 7 7^ 0. 
Without loss of generality we take 7 > hereafter. Un- 
der these conditions, the system will reach the steady 
state solution P oc (E, I) already given in ([|). The asymp- 
totic stress a is found by averaging the local strain I 
over P 0o (E,l) ) i.e. a — (I), which defines the flow curve 
(7(7). Example flow curves for different values of x are 
shown in Fig. [l| In all cases the gradient dcr/d7 decreases 
with increasing 7, indicating that the apparent viscosity 
rj = a/j decreases with 7 and that the system is every- 
where shear thinning. 

It should also be noted that there is a a one-to-one 
correspondence between 7 and a. This means that every 
point on the curve can be reached in one of two ways : 
either by applying a known strain 7(f) ~ jt, as described 
above, or by imposing a constant stress a. The unique- 
ness of the steady state solution ensures that the same 
Poo(E,l) will be reached in both cases. Of course, this 
assumes that the steady state is reached at all, for which 
it must both exist and be stable. For the standard SGR 
model, the steady state solution is stable as long as it 
exists but this does not hold for all x(a,l) under an 
imposed stress, as will be discussed in later sections. 

The behaviour of a as 7 — ► + depends upon the choice 
of x 0. For instance, for x > 2 the stress scales as a ~ 7 
and the system is Newtonian, whereas there is a power 
law fluid regime a ~ j x ^ 1 for 1 < x < 2. However, the 
most important regime for our purposes is a: < 1, for 
which er approaches a finite value. The yield stress ay, 
defined by 

cry = lim a(j) , (6) 

7-»0+ 

smoothly tends to zero as x — > 1~ and remains zero for all 
x > 1. If a system with x < 1 has a stress a < o~y applied 
to it, then it will never reach a steady state with a finite 7, 
simply because Poo(E,l) for these parameters does not 
exist. Instead the strain will logarithmically creep ac- 
cording to j(t) ~ hit ||. There may be a crossover to a 
different behaviour at late times if a is close to ay , but 
it must always be true that 7 — > + . 

The logarithmic creep in 7(4) under an imposed stress 
(7 < ay is an important result. If it is realised in 
an experimental situation, then when the strain rate 
7 ~ 1/t drops below the resolution of the apparatus, 
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which it must do at long times, it might be erroneously 
deduced that the system has stopped flowing altogether, 
i.e. 7 f=a 0. This is in accord with our intuitive notions 
about jamming : the sample initially flows when pushed, 
but later stops flowing, or 'jams.' One might argue that, 
technically, a jammed system should have 7 = exactly, 
but it is not possible for any model in this class to sus- 
tain a finite stress without a finite 7 (unless the sample 
is allowed to age for an arbitrarily long period of time 
before being sheared; see ||). This is because each indi- 
vidual element has a finite relaxation time until it yields 
and releases its stress, which can only be balanced by an 
increasing strain. 

With this insight, we now define a 'jammed' state for 
this class of systems to be a configuration which has a 
finite limiting stress ay when it is driven by an arbitrar- 
ily small strain rate 7 = + . For the SGR model, this 
corresponds to x < 1; however, for x = x(a) it may 
also depend on the history of the sample, as will now be 
discussed. 



IV. JAMMING MODEL A: x = x(a) 

It is useful to recall the physical picture that under- 
lies an x = x(a). As the system becomes sheared, either 
by an imposed strain 7 or an imposed stress a, it will 
become distorted, and this may affect its susceptibility 
to noise. The precise manner in which this happens will 
depend on the microscopic composition of the material 
in question. For the systems we are interested in, i.e. 
those that can shear-thicken or jam, the distorted con- 
figuration is less susceptible to noise that the undistorted 
one. One way to incorporate this effect is to regard the 
system as becoming 'colder' as its shear stress increases, 
which corresponds to a decreasing x{a). Thus setting 
x = x(cr) provides a mechanism for allowing the effec- 
tive temperature x to evolve with the global state of the 
system. 

In this section we shall consider x(a) that takes val- 
ues greater than 1 for some ranges of cr, and less than 
1 for others. According to the results of the previous 
section, this should allow the system to change from a 
jammed to a flowing state in response to the driving, or 
more precisely, in response to changes in a(t). Thus we 
may be able to observe a shear-induced jamming transi- 
tion (something that is not possible in the SGR model, 
for which the existence of a yield stress depends purely 
on the choice of the external parameter x). For clarity, 
we shall restrict our attention to the simplest choice of 
x(a) that has the potential to shear-thicken and jam, 
namely x(a) > 1 for small cr and x(cr) < 1 for large cr, 
with a monotonic (possibly discontinuous) behaviour for 
intermediate cr. 



A. Steady state behaviour 

Although allowing x to vary in time according to 
x = x[<j(t)] can complicate the transient behaviour of 
the system, as described below, the steady state is easy 
to analyse. This is because the very definition of a steady 
state means that a asymptotically approaches a con- 
stant value, and thus x also becomes constant. There 
is therefore a straightforward procedure for generating 
the x = x(a) flow curves from those of constant x : for 
any given value of cr, one simply finds the SGR flow curve 
for the corresponding value of x(a) and reads off the re- 
quired value of 7. This amounts to interpolating between 
the various SGR flow curves. Some examples are given 
in Fig. |^ for various x(a) that change from 1.5 for small 
a to 0.5 for large cr. 

If x(a) takes values greater than 1 for small cr, but 
smaller than 1 for larger cr, then it may jam under an 
applied stress. However, for this situation to be realised, 
the applied stress must be simultaneously large enough 
that x(a) < 1, and also small enough that a < cry • Since 
the yield stress cry also depends on x, the precise crite- 
rion for the upper yield stress to be attainable is that 
there is a range of a for which 

cr < cty[2;(ct)] . (7) 

This corresponds to the region in which the x(a) curve 
falls below the cry(x) curve when they are both plotted 
on the same axes, as in Fig. ^. For example, of the flow 
curves already presented in Fig. |], only the third example 
obeys the criterion (0) for a finite range of a. Imposing 
a stress cr from within this range will result in a system 
that is jammed in the sense that the strain logarithmi- 
cally creeps, ^(t) ~ hxt and ^{t — > 00) = 0. 

Two additional complications arise when the system is 
driven by an imposed strain rate 7 rather than an im- 
posed stress. Firstly, if 7 is increased at a sufficiently 
slow rate that the system is always arbitrarily close to 
its steady state, the stress may undergo a discontinuous 
jump from one branch of the flow curve to another, as 
marked in Figs. ||(b) and ||(c). If 7 is decreased in a sim- 
ilar fashion, then a(t) will follow the upper branch until 
it again jumps discontinuously at a different value of 7. 
Thus the system exhibits hysteresis. Hysteresis due to 
non-monotonic flow curves has also been observed exper- 
imentally ]l5|- |l7j , |l9t . The complementary form of non- 
monotonicity, which would allow a discontinuous jump 
in 7 for a small change in cr, cannot be realised in this 
class of models. This is because there is only one value of 
x(a), and hence one steady state solution, for any given 
value of cr. 

The second complication concerns the stability of the 
jammed state. In Fig. ||, the yield stress (Jy(x) and an 
x(a) that obeys (^) for a range of cr have been plotted on 
the same axes. Also plotted is the line of a that can be 
reached for a small but finite 7 > 0, which converges to 
cry (a;) as 7 — > 0. These are the a that can be realised in 
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practice, since the steady state solution (^) is not denned 
for 7 = 0. For this example there are 3 roots : a 'flowing' 
root with a = + , and two 'jammed' roots with finite 
ay > 0. The arrows in this diagram point in the direc- 
tion in which the stress will be varying for a given point 
on the line x(a). This direction is based on the assump- 
tion that, for time scales over which x can be treated as 
a constant, which will certainly include the asymptotic 
regime, the system will behave like the SGR model and a 
will evolve towards ay ■ It is clear from this that all points 
tend to flow away from the middle root, which is there- 
fore unstable. If the system is initially placed near to this 
root, it will undergo a transient and converge to either 
the higher root, or to the 'flowing root' with a = + . 

Fig. H can also be used to predict the stability of the 
steady state for finite 7. As 7 increases, the asymptotic 
stress (7(7) for the SGR model will move away from the 
yield stress curve, but will always remain a monotonic 
decreasing function of 7 [Q. Thus a smoothly decreas- 
ing x(a) will only intersect with it 1 or 3 times; when 
there are 3 roots, the middle one is always unstable, 
for precisely the same reasons as given in the previous 
paragraph. The range of a which are unstable and can- 
not be realised under an imposed 7 have been labelled 
'A' in Figs. |^(b) and (c). Note that the unstable roots 
correspond to regions of the flow curve with a negative 
gradient, as in experiments and from other theoretical 
considerations Jlj|-|ll|. There are no stability issues for 
an imposed a, for which there is always a single, stable 
root. 



B. Transient behaviour under an imposed 7 



an earlier perturbation of finite duration does not decay 
to zero with time |s|]23]-|25|| ■ 

The choice of x[a) used in the previous example is 
in fact the same as that plotted in the stability dia- 
gram, Fig. 0. This allows for a striking illustration of 
the instability of the middle root that lies near the point 
a w 0.3 on the stability diagram. Plotted in Fig. |] are 
the a(t) curves for a fixed 7 = 10~ 3 and a range of val- 
ues 0.1 < 70 < 1. It is clear that the stress will always 
diverge away from the unstable root at late times, no 
matter how finely 70 is tuned. Thus only the roots at 
a = + and a w 0.65 are stable, as previously claimed. 

There are many other complications that can arise due 
transient behaviour in a strain-controlled system. For in- 
stance, all of the a(t) curves plotted in Fig. Breach their 
global minimum cr min at a time 1 <C t <C I/7. It can be 
shown that er m i n becomes arbitrarily small as 7 — > + . 
Since this corresponds to a state with a high effective 
temperature x(<7 m i n = + ) > 1 in which the stress is 
rapidly dissipated, then if 7 is sufficiently small, a(t) will 
remain low and the system will crossover to the flowing 
root, irrespective of 70. However, this effect can be re- 
versed if the initial step shear is replaced by a smoothly 
varying j(t), such as one that exponentially decays to- 
wards its final value of 7, which is more like what could 
be attained in an actual rheometer. We will not pursue 
these complications any further here. 

In summary, the SGR model generalised to allow the 
effective temperature x to vary with the global stress a 
can exhibit, for suitable choices of x(a~) : hysteresis in 
a~(t) for slowly varying 7, as seen from the flow curves in 
Fig. ^; shear induced jamming, as in Fig. ^; and strong 
history-dependence of the existence of a yield stress. 



There are in principle no difficulties in extending the 
results described above to a system that is driven by 
a time-dependent imposed stress a(t). As long as a(t) 
tends to a constant value cfq at long times, then the same 
steady state behaviour as previously discussed will apply, 
with a replaced by do . The transient behaviour does 
not affect the final state. However, this is not the case 
when the system is instead driven by an imposed strain 
rate "f{t). In this case, a(t) can vary in a manner that 
is difficult to predict in advance, so that the final stress 
reached, and hence whether or not the system is jammed, 
will in general depend on the entire history of the system. 

A concrete example of the history dependence of a yield 
stress is given in Fig. [|. Here, the system is first subjected 
to a step shear of magnitude 70, and is then continuously 
sheared at a rate 7, i.e. jit) = jo+jt- As 7 tends to zero, 
the stress approaches an asymptotic value ay; however, 
whether or not ay is finite depends on 70. For 70 = 1 the 
stress is seen to be converging to a finite value a « 0.65, 
but for 70 = 0.1 it rather tends to a 'flowing' state with 
a = + (more precisely, a oc 7 5 as 7 — > + ). The system 
can be said to exhibit strong long-term memory, where 
the use of the word 'strong' means that the memory of 



V. JAMMING MODEL B: x = x(l) 

The second limiting form for x = x(a, I) to be consid- 
ered is one that depends only on the local strain, x = x(l). 
This marks a more drastic departure from the SGR model 
than the x(a) considered in the previous section, since 
now every element will generally have a different effec- 
tive temperature x. The steady state will therefore be 
different to that of the SGR model, and it will not now be 
possible to map the flow curves for x(l) across from those 
for constant x. One could also envisage regions of the pa- 
rameter space for which the steady state cannot even be 
reached. Indeed, this is precisely what we find : for some 
finite region of parameter space, the system reaches an 
oscillatory regime under an imposed stress, but not un- 
der an imposed strain. This is the central finding of this 
section. 

The physical justification in choosing x = x(l) is sim- 
ilar to that already discussed for x = x(a), i.e. it is 
assumed that the material can become more or less sus- 
ceptible to noise in its strained state. The main difference 
is that this is now assumed to be a local effect, that can 
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be described purely on the level of individual elements. 
Just as in the previous section, we shall focus our atten- 
tion on x(l) that decrease with I, since it is these that 
have the potential to exhibit shear thickening. 



A. Monotonicity of the flow curves 

Given x(l), the steady state stress a for a given 7 can 
be found by calculating the mean strain (I) for the steady 
state solution (||) . Some example flow curves are plotted 
in Fig. ^. They are clearly monotonic : their gradient 
is everywhere positive, and there is a one-to-one rela- 
tionship between a and 7. In these figures, the x(l) were 
chosen to vary in a stepwise manner, taking a value x > 1 
for / < 0.4 and a lower value x < 1 for I > 0.4. This is 
typical of the x(l) employed in this section. However, the 
monotonicity result is entirely general and applies for any 
x = x(Z), as demonstrated in Appendix [b|. 

The physical reason for the monotonicity of the flow 
curves is that the elements are uncoupled when x = x(l) 
and the system is driven by an imposed strain 7(t). That 
is, the expected strain reached before a given element 
yields does not depend on the state of any other ele- 
ments. This is not true for the general case x = x{o~, I), 
when the value of x for an element depends on the global 
stress cr, and thus on the state of the rest of the sys- 
tem. As long as independence holds, the average strain 
reached before each individual element yields, and hence 
the global stress a, can only increase with increasing 7. 
Thus the flow curves must be monotonic. 

The monotonicity of the flow curves is an important 
finding. As mentioned in the introduction, rheological 
instabilities often arise when the system has been driven 
to a point on its flow curve which has a negative gradient. 
Fluctuations can then cause spatial inhomogeneities to 
arise, such as shear bands with either the same stress and 
different strain rates, or equal strain rates but differing 
stresses jl^Jbifl. Alternatively, temporal oscillations may 
be observed ]l7| , p^| . Since the flow curves for x — x{l) 
have no regions of negative slope, one would not expect 
any instabilities to appear in this model. Nonetheless 
temporal oscillations in 7 do occur for a finite range of 
imposed cr, as we now describe. 



B. Oscillatory behaviour under an imposed a 

The monotonicity result described above was at- 
tributed to the independence of the elements, which is 
only true for a strain controlled system. By contrast, 
the elements become coupled under an imposed a since, 
for example, a single element yielding causes the mean 
strain to decrease slightly, which must be immediately 
countered by an increase in 7, which affects every element 
in the system. Thus collective behaviour can now occur. 
In fact, this collective behaviour cannot alter any system 



that has already reached steady flow, which we know is 
identical for strain and stress controlled systems. Thus as 
long as the stress controlled system reaches steady flow, 
it will fall onto the same monotonic flow curve as before. 
However, the couplings between the elements can dras- 
tically alter the transient behaviour, to such an extent 
that steady flow may never be reached. 

An example of collective behaviour is given in Fig. [?|, 
which shows 7(f) for a system driven by an imposed stress 
a = 0.05. As before, this plot was generated by numer- 
ically integrating the master equation ([!]) from an ini- 
tially unstrained state, using the procedure described in 
Appendix |A|. For this example, the system undergoes a 
short transient before entering into an oscillatory regime, 
in which j(t) varies with a single period of oscillation. 
There is no suggestion of a decay in the amplitude of 
the oscillation in j(t) over the largest simulation times 
attainable, even when plotted against hit (not shown), 
suggesting that this is the true asymptotic behaviour. 
For different choices of x(l) it is not possible to identify 
a single period of oscillation, as demonstrated in Fig. ||. 
It is not yet clear if this behaviour represents a distinct 
regime with more than one period of oscillation, possibly 
even a precursor to chaotic behaviour, or if it is merely a 
long-lived transient that eventually crosses over to either 
steady flow or a single-period oscillation at later times. 

The extensive simulation times means that we have so 
far been unable to map out the class of x(l) that give rise 
to an oscillatory or otherwise non-steady state. Nonethe- 
less we can make the following observation. For all of 
the oscillatory cases that we have observed thus far, x(l) 
changes from a value in the range 1 < x < 2 for small I, 
to a value x < 1 for larger I. Other choices of x(l) may 
produce an oscillatory transient, but its amplitude always 
seems to decay to zero in time. It is not clear to us why 
only this class of x(l) can produce a stable oscillatory 
regime, but it is tempting to note that the requirement 
that 1 < x < 2 for small I is also the range of x for which 
the SGR model does not have a linear regime under an 
imposed a ||. That is, a significant proportion of ele- 
ments will eventually gain a finite strain, no matter how 
small cr is, and thus the variation in x(l) will eventually 
be 'felt.' 

Once a suitable x(l) has been found, it is possible to 
move in and out of the non-steady regime by varying the 
applied stress a. Generally, we find that a high a gives 
rise to steady flow, and low a produces oscillations; how- 
ever, the excessive simulation times needed for small a 
means that we have not been able to rule out another 
crossover to steady flow as a — > + . Within the oscilla- 
tory regime, the mean strain rate averaged over a period 
of oscillation is generally much lower than that predicted 
by the flow curve. For the examples already given, the 
mean strain rate deviates from the flow curve by a fac- 
tor of 2 for the oscillations shown in Fig. |?], and by two 
orders of magnitude for the (possibly transient) oscilla- 
tions of Fig. ||. Again, we have been unable to fully 
characterise this behaviour with the available simulation 
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resources. To an extent, these problems will be alleviated 
by considering a simplified model considered in Sec. [v| . 
Before turning to consider this, we provide a qualitative 
description of the oscillatory behaviour in terms of the 
local strains I. 



C. Qualitative description of the evolution of the 
system in the oscillatory phase 



The mechanism behind the oscillatory behaviour can 
be qualitatively understood by inspecting the evolution 
of the P(E, I, t) distribution during a single period of os- 
cillation. As an example, a suitable sequence of snapshots 
is given in Fig. |{]. To aid in the interpretation of these 
figures, it is useful to recall some properties of the master 
equation (|l|) . Firstly, the convective term jdP/ dl means 
that any maxima or minima in P(E, I, t) will move in the 
direction of increasing I at a rate 7(f) . These extrema will 
eventually disappear when large values of I are reached, 
and the yield rate T(E, I) becomes very high. Also, there 
is a flux of newly-yielded elements at / = 0, which have 
barriers E weighted according to p(E). 

Since the oscillatory behaviour is by definition cyclical, 
it is somewhat arbitrary where one chooses to start the 
sequence of snapshots. In Fig. || we have chosen times 
corresponding to before, during and just after the point 
when 7 reaches its maximum value. In the first snapshot, 
P(E, I, t) is concentrated into two regions : a sharp peak 
at I ~ 0, and a broader peak with I in the range 1 < I < 3. 
The hrst peak is 'hot' in the sense that it has the higher 
value of x(l), but since it has a low strain, it does not 
significantly contribute to the total stress a. The second 
peak, although highly strained, lies in a region in which 
x(l) is small, and therefore has a low rate of yielding. As 
long as the yield rate is low, so too is the rate of stress 
loss from elements that belong to this peak. Therefore 
the strain rate 7 will also be low, and indeed this first 
example corresponds to a system in which 7 is close to 
its minimum value. 

This state of affairs is not stable, however. Although 
the yield rate of the highly strained elements is low, it 
is nonetheless non-zero, and therefore so is 7. Conse- 
quently the whole system is being convected at a finite 
rate, so the elements in the high-Z region are becoming 
more strained. This decreases their effective energy bar- 
rier AE = E — \l 2 , which increases their yield rate. But 
an increased yield rate means an increased 7, which in 
turn increases the rate at which the elements are becom- 
ing more strained, which increases their yield rate yet 
more, and so on. This description is that of a positive 
feedback loop, which causes 7 to increase at ever faster 
rates until all of the heavily strained elements have been 
depleted. The second snapshot in Fig. ^| show the state 
of the system shortly before this has happened. 

At the same time that the highly-strained elements are 
being depleted, 7 is close to its maximum value and con- 



sequently P(E,l,t) in the small I region becomes some- 
what flat, certainly much flatter than under a small 7. 
This can clearly be seen in the second and third snap- 
shots of Fig. ^. When 7 again falls to lower values, this 
flat part of the distribution will start to decay as the el- 
ements within it yield. However, the yield rate depends 
on x(l), and since x(l) changes to a lower value at I = 0.4, 
P(E, I, t) will decay much more rapidly for I just smaller 
than 0.4 than for I just greater than 0.4. Thus a dip will 
occur around the point I ~ 0.4, which can clearly be seen 
in the final snapshot. As time increases, this dip will 
become more pronounced until the system can again be 
separated into a population of highly-strained elements 
and a second population with I w 0. This is where we 
began, and thus the cycle is complete. 



VI. JAMMING MODEL C: x = x(l), 
MONODISPERSE E 

It is possible to more completely describe the oscilla- 
tory regime for a simplified version of the model where 
every element has the same energy barrier E\ . This 
formally corresponds to a 'monodisperse' prior distribu- 
tion p(E) = 5{E — Ei), as opposed to the exponential 
p{E) = c~ E which has been employed in all of the ear- 
lier sections. The reduction in the number of degrees of 
freedom that this entails significantly decreases the sim- 
ulation times, and therefore allows for a more thorough 
numerical investigation of the model. It is also possible 
to derive an analytical criterion for the stability of the 
steady state, as described below. 

Furthermore, the very fact that a monodisperse system 
also has an oscillatory regime clearly demonstrates that 
this phenomenon is robust. In particular, it does not re- 
quire the usual SGR assumption of an exponential tail to 
p(E). Therefore the possible existence of a yield stress, 
which was so central to the history-dependent jamming 
scenario at controlled 7 described in Sec. IV, is not nec- 
essary. This robustness means that the mechanism be- 
hind the oscillations at controlled a may be realisable 
in a much broader range of models than the restricted 
set considered here and, possibly, may also occur in real 
materials. 

The emergence of oscillations in this monodisperse 
model has been separately discussed elsewhere |2(| . Here 
we supply extra details, such as the stability analysis of 
the steady flow state. For the sake of completeness, the 
overall behaviour of the model will also be briefly de- 
scribed, although we refer the reader to [^0| for a fuller 
discussion of many of the points raised below. 



A. Steady state behaviour for monodisperse E 

Perhaps the most immediate consequence of only al- 
lowing a single barrier E\ is that the system can never 



8 



have a finite yield stress. Indeed, in the linear regime 
7 — > + , the steady state stress is just a ~ -ye 5 ' 1 /^ ), 
which vanishes with 7. There are no qualitative differ- 
ences for ar<l,l<x<2 etc., as in the SCR model. 
The complete lack of a yield stress means that a monodis- 
perse system with x — x(a) will not exhibit a jamming 
transition, in contrast to the situation for a p{E) with an 
exponential tail described Sec. IV. 

One respect in which the monodisperse E model is sim- 
ilar to the exponential p(E) case is that its flow curves 
for x — x(l) are also monotonic. In fact, the monotonic- 
ity proof given in Appendix [b] holds for arbitrary x(l) 
and p(E). Thus there are no regions on the flow curve 
with a negative slope, and therefore there are no ranges 
of control parameters for which one might normally ex- 
pect bulk shear banding to arise. Just as in the previous 
section, however, oscillatory behaviour can be realised 
under an imposed constant stress. 



B. The oscillatory regime 

As in the polydisperse case, the monodisperse model 
exhibits an oscillatory regime under conditions of im- 
posed stress, but not under an imposed 7. However, the 
oscillatory regime in the monodisperse model differs from 
the polydisperse case in that only single-period oscilla- 
tions have so far been observed. There is no suggestion of 
the more complex non-steady behaviour hinted at earlier 
in Fig. ||, for example. Some examples of the oscillatory 
behaviour in the monodisperse case are given in Fig. |lC| . 
Here, the strain -f(t) is shown for 3 systems with Ex = 5 
and the same x(l), but different imposed a. Remarkably, 
the mean strain rate (7), where the use of the angled 
brackets now denotes 7 averaged over a single period of 
oscillation, is clearly a decreasing function of a, in com- 
plete contrast to the monotonic flow curve. 

Varying a over a wider range of values reveals that 
steady flow is reached for either sufficiently small or suf- 
ficiently large imposed stresses; only for intermediate a 
are oscillations observed. The (7) as read off from the 
simulations are plotted against a in Fig. O, overlayed 
with the steady state flow curve. It is clear that, if the 
system reaches a steady state, it falls onto the flow curve 
and therefore the dependence of (7) with a is monotonic. 
Within the oscillatory regime, however, the (7) line devi- 
ates from the flow curve to an extent that it even becomes 
non-monotonic for a broad range of a. 

The shape of the oscillations varies with a. Close to 
either transition to the steady flow regime, the oscilla- 
tions are approximately sinusoidal, indicating that there 
is only a single, finite frequency of oscillation at the tran- 
sition, the amplitude of which vanishes with the onset of 
steady flow. Further into the oscillatory regime, 7 can no 
longer be decomposed into a single harmonic, but instead 
approaches a waveform in which most of the variation in 
7 is compressed into a small fraction of the total period 



of oscillation. Consequently j(t) approaches an almost 
rectangular, staircase shape. The mechanism underly- 
ing the rapid increase in 7 is the same positive feedback 
loop that has already been discussed for the polydisperse 
case, as can be seen by ins pec ting the evolution of the 
P(l) distribution with time |20[ . 

Near to the transition points at a = o~ c , the ampli- 
tude of the oscillation appears to vanish as \a — a c \ a 
with a > 0, as demonstrated in Fig. [jj]. Data fitting 
suggests that the lower threshold lies at a c w 0.07 with 
a value of a w 0.7, and that the upper threshold is at 
o~ c w 1.204 and has a different value of a rj 0.2. How- 
ever, the range of values given in this figure is somewhat 
limited, constituting only an order of magnitude of vari- 
ation in \o~ — <j c | and an even narrower range of ampli- 
tudes. The reasons for this are purely technical : close to 
the transition points, the amplitude decays very slowly 
in time to its asymptotic value, rendering the required 
simulation times impractically long. Hence it is conceiv- 
able that the true values of a are significantly different 
from those found here. In particular, a value of a = 0.5 
for both cases, as expected for a Hopf bifurcation [ p6[ , 
cannot be ruled out. 

Finally, plotted in Fig. [l3] is the product of the mean 
strain rate (7} and the period of oscillation T for different 
values of a. To first order, {y)T is seen to be independent 
of (j, which suggests that the oscillatory regime can be 
characterised by a single strain I* ~ (j)T. In this case 
I* w 2.3, which is also the typical increase in 7 during 
a single cycle of oscillation. A possible interpretation of 
I* is that it is the amount by which the system needs 
to be strained until the positive feedback loop discussed 
earlier starts to dominate the system behaviour, causing 
it to 'reset' to the start of its cycle. If this is correct, 
then I* should correspond to the point at which highly- 
strained elements have the same yield rate as unstrained 
elements, i.e. T(0) ~ T(l*). Rough calculations based on 
this assumption give I* ~ ^/2£'i[l — x(oo)/x(0)], which 
for this example predicts I* « 2.4, in fair agreement with 
the observed value. 



C. Analysis of the stability of the steady state 

A second advantage of the monodisperse model is that 
it is possible to derive an analytical criterion for the sta- 
bility of the steady state at fixed driving shear stress a, 
although even in this simplified case the resulting expres- 
sion is still somewhat unwieldy. Suppose steady flow is 
reached with a strain rate 700 = 700 (ct) . Then the corre- 
sponding Poo (I) is 



Poo(0 =AAexp 



-7— — fv{l')Al' 

7oo Jo 



(8) 



where the normalisation constant J\f is fixed by the con- 
straint / °° Poo (I) dl — 1. We now look for eigenfunctions 
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{p(l),g} and eigenvalues {s} of the linearised relaxation 
operator by writing 



P(/,i) = P oo (0 + s P (l)e s 



7(t) 



7o 



(9) 
(10) 



In these expressions, g is a real constant, and the real 
functions p(l) remain to be found. The complex con- 
stant s determines the stability of the steady state : it 
is unstable precisely when Re(s) > 0, since the ampli- 
tude of the perturbation will then increase exponentially 
in time. Conversely, it is stable when Re(s) < 0. The 
frequency of the oscillatory part of the motion near to 
the steady state is Im(s)/27r |2(|. 

The unknown function p(l) can be found by insert- 
ing (H) and ([h]) into the master equation and neglecting 
terms of 0{t^), which gives 



p'{l)i 00 +p{l)[s + T(l)]=-gP' 00 {l) 

_ <?r(/) 



Too 



Poo(0 



(11) 

(12) 



This is just a first-order differential equation in I and 
integrates to 



p(l) = — 

Too 



A 



x exp 



- / d; 

o 



dl'T(l')e sl '^- 
fS + T(l>) 



(13) 



The constant A can be found by imposing dl P(l, t) — 
1, i.e. f™dlp(l) = 0. 

Also, since we are considering a stress-controlled sys- 
tem, a = (I) must remain constant and hence 



dlp(l)l = 



(14) 



This global constraint allows the value of s to be speci- 
fied. Inserting (]l|) into ^ gives, after some manipula- 
tions, the following equation for s, 



oo poo 



dh / dh (h - h)PUh)Poo{h) e- s ^+ 1 ^/^ 
o Jo 

x f 2 dir(l)e sl/ ^°° =0 (15) 
Jo 

Although this equation cannot be rearranged to give an 
expression for s in terms of the system parameters, as 
would be desired, it is still enough to show that a purely 
exponential variation away from the steady flow cannot 
be observed. This follows from the observation that the 
first line in ( |l5| ) is an odd function in 1% — li, but the 
second line is a strictly increasing function of li when 
Im(s) = 0. Hence the equation cannot be obeyed for 
a purely real s, and the transient behaviour close the 
the steady state must include an oscillatory component, 
which is consistent with the simulation results. 



VII. DISCUSSION 



One of the most striking findings of this work has been 
the observation of oscillations in j(t) under a constant 
imposed stress a for some x — x(l), as described in 
Sees. [v| and VI. Some consequences of this phenomenon 
have already been discussed in p0| |. Two further points 
will be discussed here : the identification of the dominant 
physical mechanisms, and the relationship to so-called 
'stick-slip' behaviour observed in other systems. Both of 
these will be considered in turn. 



A. Physical picture and analogous phenomena 



In Sec. V C , the mechanism behind the oscillatory be- 
haviour was explained in terms of two separated popula- 
tions of elements : a 'cold' population of highly strained 
elements with a low x, and a 'hot' population of elements 
with / w and a high x. It was explained how the cold 
elements can give rise to a positive feedback loop, caus- 
ing j(t) to rapidly increase until the cold elements have 
yielded. At the same time, a new population of cold ele- 
ments is produced from the hot one. 

Known instances of rheological instability are often ex- 
plained in terms of the spatial coexistence of subpopu- 
lations, or phases. For instance, the temporal oscilla- 
tions in viscosity observed in wormlike micelles under 
an imposed stress was attributed to a slowly fluctuat- 
ing interface between a fluid phase and shear-induced 
structures [fill. For surfactant solutions in the lyotropic 
lamellar phase, it was attributed to coexisting ordered 
and disordered phases ||l7| . However, these instabilities 
occur in the vicinity of non-monotonic regions of the flow 
curve. Spatial instabilities, such as shear banding, can 
also be found near to where the flow curve has a nega- 
tive gradient By contrast, the temporal oscillations 
observed in our models arise even though they are by as- 
sumption spatially homogeneous, and the flow curves are 
everywhere monotonic. 

This suggests that the mechanism behind the oscilla- 
tory behaviour seen here has not yet been observed in a 
rheological context. It is therefore sensible to look further 
afield for analogous phenomena. One possibility comes 
from mathematical biology. The Glass-Mackey equation 
describes the variation in time of the size of a population 
of white blood cells in response to a hormonal control 
system [£7| . It has a form similar to that of a first-order 
differential equation, but includes a feedback term that 
depends on the state variable at an earlier time. This 
delay corresponds to the maturation time of white blood 
cells from stem cells. Because of this delay, the pop- 
ulation size can vary in time in a non-trivial manner, 
including oscillatory behaviour and chaos M . 
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That some form of time delay may also play a role 
within our x(l) model is clear : once an element becomes 
'cold,' its yield rate remains very low until some time 
later, when it becomes sufficiently strained that its yield 
rate becomes comparable to that of 'hot' elements. Thus 
there is a delay between when an element first becomes 
cold and when it yields, although here the delay time 
is not constant but depends on j(t). Thus it is possi- 
ble that the oscillatory behaviour observed for x = x{l) 
could be described by a simplified equation, similar in 
form to the Glass-Mackey equation. This is a particu- 
larly enticing proposition, as the Glass-Mackey equation 
is capable of producing chaos, and chaotic behaviour has 
also been observed in surfactant solutions p8 29|, albeit 
in the signature of o~(t) under an imposed 7. It is not 
yet clear to us if a meaningful mapping between the two 
models is possible. 



B. Oscillatory behaviour or 'stick— slip' ? 

Deep into the oscillatory regime, the waveform of j(t) 
throughout a single period of oscillation becomes increas- 
ingly 'sharp,' with most of the variation in j(t) occurring 
in just a small fraction of the total period of oscillation. 
The ratio of the maximum value of j(t) to its minimum 
has been shown to exceed two orders of magnitude, and 
we see no reason to deny that greater separations may 
be attained for different parameter values. 

It is tempting to refer to this behaviour of j(t) as 
'stick-slip' behaviour, in which the system is 'stuck' until 
the short duration of time in which 7 rapidly accelerates 
to its maximum value, corresponding to a 'slip' event 
This kind of response is also observed in earthquakes 
ultra-thin liquid films |l(J and granular media |]ll| , |l2 
for example. However, we have instead chosen to refer to 
the variation in j(t) as merely 'oscillatory,' since the un- 
derlying physics seems to be somewhat different to the 
examples of stick-slip behaviour mentioned above. In 
particular, the term stick-slip is usually employed to re- 
fer to a surface friction phenomenon, whereas the models 
studied in this paper have no surface, or indeed any form 
of spatial definition. They are only intended to describe 
the bulk behaviour of a material. Furthermore, we only 
find oscillations in 7(f) under an imposed er, and never 
vice versa, whereas stick-slip seems to more usually (al- 
though not exclusively) refer to variations in the (nor- 
mal) stress under a constant driving velocity. Thus to 
avoid possible confusion, we choose not to refer to the 
behaviour observed in the models as stick-slip. 



of Sollich et ai, but differ in that the effective tempera- 
ture x is no longer constant, and can instead vary with 
the state of the system through either the global stress 
a or the local strain I. We have considered choices of 
x that decrease for increasing er or I, which is relevant 
to the study of shear thickening materials. These mod- 
els have no spatial definition, and thus by construction 
cannot exhibit any form of spatial heterogeneity. 

For x = x(a), the flow curves can be extracted from 
the known curves for constant x. Many choices of x(o~) 
produce flow curves with non-monotonic regions, which 
exhibit hysteresis in o~(t) under ramping the strain rate 
j(t) first upwards and then downwards. Furthermore, 
a subset of these x(cr) also give rise to a jammed state 
for a range of applied stresses, in that the strain 7(4) 
creeps logarithmically, j(t) oc ln(t). The criterion for 
this to arise is that the curve of x(a) drops below the 
SGR yield stress curve try (%) when they are plotted on 
the same axes. For an imposed strain rate that decays to 
zero at late times, a jammed configuration was defined 
as one with a finite asymptotic stress, a(t) ~ ay > as 
j(t) — > + and t — > 00. It was found that whether or 
not a jammed configuration was reached depends on the 
entire strain history of the system, a situation that was 
referred to as history dependent jamming. 

For x = x(l), the flow curves are always monotonic, 
and steady flow is always reached under a constant im- 
posed strain rate j(t) ~ jt. However, for a range of im- 
posed stresses and some choices of x(l), steady flow is not 
realised. Numerical integration of the master equation 
demonstrated that j(t) oscillated around a well-defined 
mean with a single period of oscillation. The possibility 
of more complex non-steady behaviour in some regions 
of parameter space could not be ruled out. A similar os- 
cillatory behaviour occurs with a simpler model in which 
every element has the same energy barrier |20 , which 



suggests that this phenomenon is robust. Finally, we dis- 
cussed the relationship between this oscillatory behaviour 
and that observed in experiments, and considered analo- 
gous phenomena from fields outside of rheology. 
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VIII. SUMMARY 

We have introduced a range of schematic models that 
are capable of exhibiting a form of jamming under an im- 
posed stress a. The models are based on the SGR model 



APPENDIX A: SIMULATION DETAILS 

Direct numerical integration of the master equation (|l|) 
has proven to be unstable with respect to discretisation 
errors. Instead, the results in this paper were generated 
from the transformed equation 
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dP{E,Al,t) 
dt 



= -T(Al + 1 (t))P(E,Al,t) 
+ oj(t) 5(Al + j(t)) p(E) 



(Al) 



where = I — 7(t). This removes the convective term 
and dramatically improves numerical stability. The dis- 
crete probability distribution = P(iSE, jSl) was de- 
fined on a rectangular mesh of points {ij} and evolved ac- 
cording to ( Al) by using a straightforward Euler method. 
A further refinement was to average both the evolution 
equation (Al) and the initial conditions over the ranges 
(E, E + 5Ef and (Al, Al + 51). This significantly reduces 
the number of mesh points required for the simulations 
to properly converge, without unduly increasing the al- 
gorithmic complexity. The delta function was treated as 
a triangle of base width 261 and height 1/51, but any 
sufficiently narrow function gave the same results. 

Two classes of initial conditions were employed, but 
the long-time behaviour of the system was found to be 
identical in both cases. Only the short term behaviour 
varied, and then only in a non-essential manner. For 
the sake of completeness, the initial conditions were : 
(i) A 'quench' Pq(E,1) = 5(l)e~ E , which corresponds to 
the unstrained equilibrium at x = oo, or (ii) Pq(E,1) = 
5(l)^e- E / E ° with£ = [l-l/x(l = 0, cr = 0)]-\ which 
corresponds to an unstrained system that has been al- 
lowed to reached equilibrium. Note that this second ini- 
tial state is only defined for x(l = 0, er = 0) > 1 p3[ . 

The strain 7(f) is only known a priori for a strain- 
controlled system. When it is rather a known stress a(t) 
that is applied, 7(f) must be chosen at every time step so 
that the actual stress does not differ from a(t) by more 
than a tolerance parameter e <C 1. To ensure that this 
condition was satisfied in our simulations, a series of trial 
strain rates 7W, 7W, <-y( 3 ) . . . were generated and tested 
on a replica mesh P*a . The Pij were not updated until a 
suitable 7 had been found. 

The trial values 7™ were generated as follows. For a 
continuous time variable, j(t) = (lT(l))p^ + a(t) ex- 
actly, as seen by multiplying the master equation (|l|) 
by I and integrating over all E and I. This is there- 
fore the sensible choice for 7W. However, integrating 
over a finite time step 5t inevitably introduces errors of 
0(5t), and so the integrated stress will differ from the 
required value by some small amount 5a. Reintegrating 
with 7^) = 7W — 5a will therefore reduce the error to a 
smaller amount 0(5t 2 ). This procedure can be repeated 
to generate a series of successively better estimates j( 3 \ 
7W etc. For our choice of e = 1CP 6 , we have found that 
typically 3 — 5 such trials are needed at every time step. 



APPENDIX B: MONOTONICITY OF THE FLOW 
CURVES 

The purpose of this appendix is to show that the flow 
curves are monotonic for any yield rate T(l) that depends 
only on the local strain I. This includes the thermally 



activated T(l) with x = x(l), possibly constant, but not 
when x also depends on the stress a. It also applies for 
an arbitrary prior barrier distribution p(E). 

The steady state solution P QO (E,l) has already been 
given in (0). The asymptotic yield rate Woo = 
lim^oo uj(t) is fixed by ensuring that this expression nor- 
malises to unity, 



dE / dl p(E) e 



(Bl) 



where we have introduced the short-hand notation 
f(a,l) = £T(a,l')dl'. The stress is a = (l)oo, where the 
angled brackets '()oo' denote the average over P 00 (E,l). 
If r depends on a, then a a must be chosen that is con- 
sistent with ([|). In general there can be more than one 
suitable a. 

First consider the case T = T(l), so / = /(£). Then by 
differentiating (l)^ and using (0), 



da 



UJ 



57 T + T^^ 



(B2) 



Similarly, the equation 
give 



can be differentiated to 



-^ = i-4</(0) oo 

Woo 07 7 7 
Combining these two expressions produces 
. 9 da 



7 



(Z/(0)oo - ff(/(0)oo 

= ((!-*) [f(l)-f(a)]) c 



(B3) 



(B4) 

(B5) 
(B6) 



The final line in this equation is valid as f(a) is a con- 
stant, so ({I - a)f(a)) 00 = f(a){l - a)^ = 0. Since f(l) 
is an increas ing function of /, the quantity in the angled 
brackets in (B6) is positive both for I < a and for I > a, 
and vanishes smoothly at I = a. Thus the left hand 
side must be positive, i.e. a always increases with 7, as 
claimed. 

The same conclusion does not hold when T = T(a, I). 
Since a depends on 7, differentiating the steady state so- 
lution now gives rise to an extra term on the right hand 
side of @, 



da 



duj a 



7 7 



-i (I g{a, I)) oa I? 
7 oj 



(B7) 



using g(a,l) — J Q [dT(a,l')/da]dV. Proceeding as before 
da ((l~a)f(a,l)) 00 



7 2 — 

1 <9 7 



l+-(lg(a,l)) c 

7 



(B8) 
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Thus it is now possible for the gradient to diverge if 
(lg(a,l)} is sufficiently negative, i.e. if r(cr, /) decreases 
sufficiently rapidly with er. For the particular case of 
r(<j, I) — e~( B ~2 l )/ x ( <T ) j the gradient of the flow curve 
diverges at any point such that 



x'(a) 



-x 2 (a) 7 



lg(E- \lv)e-( E -¥' 2 )/* dl' 



(B9) 



We can see no obvious physical interpretation of this 
mathematical criterion. 



[24] 
[25] 
[26] 

[27] 
[28] 
[291 



M. M'ezard, in Spin glasses and random fields (World 
Scientific, Singapore, 1998). 

L. F. Cugliandolo and J. Kurchan, Phys. Rev. Lett. 71, 
173 (1993). 

L. F. Cugliandolo and J. Kurchan, Phil. Mag. 71, 501 
(1995). 

P. Glendinning, Stability, instability, and chaos: an in- 
troduction to the theory of nonlinear differential equa- 
tions (Cambridge University Press, Cambridge, 1994). 
L. Glass and M.C. Mackey, From Clocks to Chaos 
(Princeton University Press, Princeton, 1988). 
R. Bandyopadhyay, G. Basappa and A.K. Sood, Phys. 
Rev. Lett. 84, 2022 (2000) 



cond- 



R. Bandyopa dhyay and A.K. Sood, preprint 
mat/00 1248^ 

C. Monthus and J.-P. Bouchaud, J. Phys. A 29, 3847 
(1996). 



M.D. Ediger, C.A. Angell and S.R. Nagel, J. Phys. C 100, 
13200 (1996). 

S. F. Edwards and D. V. Grinev, in Jamming and Rheol- 
ogy, eds. A. J. Liu an d S. R. Nagel (Taylo r and Francis, 
New York, 2001); also \cond-mat/9905114. 



A. J. Liu and S. R. Nagel, Nature 396, 21 (1998). 

V. Trappe, V. Prasad, L. Cipelletti, P. N. Segre and 

D. A. Weitz, Nature 411, 772 (2001). 

M.E. Cates, J. P. Wittmer, J.-P. Bouchaud and P. 

Claudin, Phys. Rev. Lett. 81, 1841 (1998). 

P. Sollich, F. Lequeux, P. Hebraud and M. E. Cates, 

Phys. Rev. Lett. 78, 2020 (1997). 

P. Sollich, Phys. Rev. E 58, 738 (1998). 

S. M. Fielding, P. Sollich and M. E. Cates, J. Rheol. 44, 

323 (2000). 

C.H. Holz, Nature 391, 37 (1998). 

I.S. Aranson, L.S. Tsimring and V.M. Vinokur, preprint 
cond-mat/010118b\ 

M. Lubert and A. de Ryck, Phys. Rev. E 63, 021502 
(2001). 

I. Albert, P. 
L.Barabasi 



Tegzes, R. Albert, J.G. Sample , A. 
T. Vicsek and P. SchifTer, preprint 



mat/ 0102521 



cond- 



P. D. Olmsted, Curr. Opin. Coll. Int. Sci. 4, 95 (1999). 

C. -Y. D. Lu, P. D. Olmsted and R. C. Ball, Phys. Rev. 
Lett. 84, 642 (2000). 

O. Volkova, S. Cutillas and G. Bossis, Phys. Rev. 
Lett. 82, 233 (1999). 

D. Bonn, J. Meunier, O. Greffier, A. Al-Kahwaji and 
H. Kellay, Phys. Rev. E 58, 2115 (1998). 

A. S. Wunenburger, A. Colin, J. Leng, A. Arneodo and 
D. Roux, Phys. Rev. Lett. 86, 1374 (2001). 
Y. T. Hu, P. Boltenhagen and D. J. Pine, J. Rheol. 42, 
1185 (1998). 

H.M. Laun, J. Non-Newt. Fluid Mech. 54, 87 (1994) 
D.A. Head, A . Ajdari and M.E. Cates, preprint 



mat/0108224 



cond- 



F. Lequeux and A. Ajdari, Phys. Rev. E 63, 030502 
(2001). 

J.-P. Bouchaud and M. M'ezard, J. Phys. A 30, 7997 
(1997). 

J.-P. Bouchaud, L. F. Cugliandolo, J. Kurchan and 




FIG. 1. Flow curves for a system with a constant x, i.e. 
the SGR model. From top to bottom, each line corresponds 
to a value of x increasing from 0.25 to 2.5 in steps of 0.25. 
The lines x — 1 and x = 2 have been highlighted. A finite 
yield stress <ty = lim-y^o(f) exists only for x < 1. The prior 
distribution is p(E) — e~ E . 
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FIG. 2. Flow curves for x(a) = 1 - 0.5 tanh[a(cr - 0.4)] 
(thick lines), where a = 3.5 (a), 4.5 (b) and 10 (c). For com- 
parison, the thin lines are the constant x curves from Fig. 
In (b) and (c), the vertical dashed line represents the discon- 
tinuous jump in stress for a slowly increasing 7, and the region 
marked 'A' denotes the range of a that is unstable under an 
imposed 7, but is stable under an imposed a. In (c), the 
range of a for which the system is 'jammed,' i.e. cr(0 + ) > 0, 
has also been marked. The x(a) for different a are plotted 
in (d), where it can clearly be seen that only for the a = 10 
case does x(a) drop below the yield stress curve. 
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FIG. 3. Plot of the yield stress oy{x) from the SGR model, 
overlayed with a particular choice of x(a) that changes value 
from 1.5 to 0.5 with increasing stress (the actual functional 
form is x(a) = 1 — i tanh[25(<r — 0.25)]). The dashed line is a 
schematic representation of the asymptotic stress ait — » oo) 
for a small but finite 7, which is what is actually attainable 
for these models. The arrows represent the direction in which 
the stress will vary for a constant x, and explain the given 
stability of the roots. 




10 10 10 

Y(t) - Y 

FIG. 4. Plot of stress versus strain for 7(t) = 70 + jt and 
the same x(a) as in Fig. ^. The upper set of lines correspond 
to 70 = 1 and the lower set to 70 = 0.1. Within each set, 
the lines refer to (from top to bottom) 7 = 3 x 10 -3 , 10 -3 , 



3 x 10~ 4 and 10" 



As 7 



, the upper curves are seen 
to be converging to a — > <ty ~ 0.65, whereas the lower curves 
are approaching a zero-stress state according to a ~ kj ' 5 , 
where k is an arbitrary constant. In all cases the system was 
first allowed to reach equilibrium under zero shear before the 
step shear 70 was applied. 
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FIG. 5. The variation in stress for the system of Fig 
driven by the imposed strain y(t) — 70 + yt with 7 = 
demonstrating the inability to reach the root at a w 0.3 
bottom to top on the left hand side, the 'initial condition' 70 
takes the values 70 = 0.1, 0.4, 0.48, 0.49, 0.491, 0.5, 0.6 and 1. 
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x(l) = 1.5 (l<0.4), else 0.5 
x(l) = 2.5 (l<0.4), else 0.5 
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FIG. 6. Examples of typical flow curves when x depends on 
the local strain I. Here, x(l) = 0.5 for I > 0.4, but takes the 
higher value of 1.5 (solid line) or 2.5 (dashed line) for I < 0.4. 
These lines were generated from the steady state solution of 
the master equation. 
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FIG. 7. The strain ^(t) under an imposed stress ao = 0.05 
for a system in which x depends on the local strain according 
to x(l) = 1.8 for I < 0.4 and x(l) = 0.8 for I > 0.4. (Inset) The 
strain rate 7(f) for the same data after the transient. 
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FIG. 8. The strain 7(t) under an imposed stress ao = 0.1 
for x(l) = 1.5 for I < 0.7, 0.5 for / > 0.7. Despite the long 
times attained, there is no clear indication of a single period 
of oscillation. 
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FIG. 9. Snapshots of P(E, I, i) for the same system as in 
Fig. at 3 different times. For clarity only 3 values of E are 
shown, namely E — 8 (solid line), E = 10 (dashed line) and 
E — 12 (dot-dashed line). The chosen times correspond to 
just before (a), during (b) and just after (c) the point at which 
7 reaches it maximum value. The period of the oscillation is 
approximately At = 1.6 x 10 4 , so that the time between (a) 
and (c) comprises roughly | of a single oscillation. In each 
case, the / at which x(l) changes from 1.8 to 0.8 is represented 
by a vertical dotted line, and the arrow points to the same dip 
in the E = 10 distribution, which moves to the right under 
the action of 7. 




FIG. 10. The strain j(t) ior a^monodisperse system with 
Ei = 5 under an imposed step stress cr = 0.1 (solid line), 
0.13 (dashed line) and 0.2 (dot-dashed line) at a time t = 0, 
demonstrating a decrease in the mean 7 with a. Here, x(l) = 1 
for / < 0.4 and x = 0.4 for I > 0.4. The system was initially 
unstrained. 



FIG. 11. The flow curve for the same system as in Fig. [l^ 
under an imposed stress. The circles are the 7 as measured 
from the simulations, either in the steady state (solid circles) 
or in the oscillatory regime (open circles). In the latter cases, 
7 was averaged over a whole number of oscillations. The sizes 
of the circles are larger than the error bars. For compari- 
son, the theoretical flow curve for the steady state solution 
Poc(E, I) is plotted as a solid line. 
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FIG. 12. The amplitude of the oscillation, defined as 
^|7max — 7min|, against the distance from the transition be- 
tween the steady and oscillatory regimes. The upper data set 
(circles) corresponds to the transition at a c = 1.204, and the 
lower set (triangles) corresponds to <r c = 0.07. For compari- 
son, the solid straight line has a slope of 0.21, and the dashed 
line has a slope of 0.7. 
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FIG. 13. The product of the mean strain rate (7) and the 
period of oscillation T for different stresses a, demonstrat- 
ing that this quantity is approximately constant. The solid 
horizontal line represents {j)T = 2.33. 
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